---
title: "Burkholder, Hansen and Wager Replication Code"
output:
  html_document: default
  pdf_document: default
date: "2025-05-30"
---


```{r, results='hide'}
library(cowplot)
library(ggplot2)
library(readr)
library(foreign)
library(interactions)
library(effects)
library(lmtest)
library(sandwich)
library(ggeffects)
library(stargazer)
library(scales)
library(dplyr)
library(interplot)
library(modelsummary)

# Clear the workspace
rm(list = ls())

#Setting workspace and loading cleaned data
setwd("[filepath]")

df <- read_rds("Burkholder_et_al_Replication_Data.RDS")

# Recode 
df$race_5 <- factor(df$race_5, levels = c("Other/Mixed", "White", "Black", "Hispanic", "Asian"))

#Getting rid of "Other"
df <- df %>%
  filter(race_5 != "Other/Mixed")

df$race_5 <- recode(df$race_5, "Hispanic" = "Latino")


```

```{r}
# Calculate the proportions of each race category
 df %>%
  group_by(race_5) %>%
  summarise(n = n(),
            proportion = n() / nrow(df))
```


#Computing Rank Example
```{r}
educ_rank_by_state <- df %>%
  filter(educ_c == "12", year == "2020", generation == "Baby Boomers", state_abbr == "WV") 

head(educ_rank_by_state$educ_rank_generation)

educ_rank_by_state <- df %>%
  filter(educ_c == "12", year == "2020", generation == "Baby Boomers", state_abbr == "WA") 
head(educ_rank_by_state$educ_rank_generation)

educ_rank_by_state <- df %>%
  filter(educ_c == "16", year == "2020", generation == "Millennials", state_abbr == "DC") 
head(educ_rank_by_state$educ_rank_generation)

educ_rank_by_state <- df %>%
  filter(educ_c == "16", year == "2020", generation == "Millennials", state_abbr == "MS") 
head(educ_rank_by_state$educ_rank_generation)
```



# Fig 2: Average Educational Rank Table

```{r}

df$year2 <- as.character(df$year)
df$year2 <- as.Date(paste0(df$year2, "-01-01"), "%Y-%m-%d")
df$educ_rank_decile <- cut(df$educ_rank_generation, breaks = seq(0, 1, by = 0.05), labels = FALSE, include.lowest = TRUE)

# Summarize data based on educational rank deciles
summary_data_deciles <- df %>%
  group_by(educ_rank_decile, race_5) %>%
  summarise(
    avg_educ_rank = mean(educ_rank, na.rm = TRUE),
    avg_turnout = 100 * mean(voted_c, na.rm = TRUE) # Calculate average voter turnout (%)
  ) %>%
  ungroup()  # Remove grouping for subsequent operations
# Custom theme definition

# Plot creation
p1_deciles <- ggplot(summary_data_deciles, aes(x = avg_educ_rank, y = avg_turnout, color = race_5)) +
  geom_line(size = 2) +  # Increase line thickness
  geom_point(size = 3) +   # Increase point size
  labs(x = "Average Educational Rank", y = "Average Voter Turnout (%)", color = "Race", caption = "Data: CPS Voter Supplement") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  my_theme(base_size = 14)

p1_deciles

ggsave("Educational Rank and Turnout.png", plot = p1_deciles, width = 11, height = 7)

```

# Fig 2: Average Education Year Table


```{r}

summary_data <- df %>%
  group_by(educ_c, race_5) %>%
  summarise(
    avg_educ_c = mean(educ_c, na.rm = TRUE),
    avg_turnout = 100 * mean(voted_c, na.rm = TRUE) # Calculate average voter turnout (%)
  ) %>%
  ungroup()  # Remove grouping for subsequent operations


p2 <- ggplot(summary_data, aes(x = avg_educ_c, y = avg_turnout, color = race_5)) +
  geom_line(size = 2) +  # Increase line thickness
  geom_point(size = 3) +   # Increase point size
  labs(x = "Average Educational Attainment in Years", y = "Average Voter Turnout (%)", color = "Race", caption = "Data: CPS Voter Supplement") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 100)) +
  my_theme(base_size = 14)

p2
ggsave("Average Self-Reported Turnout by Education and Race.png", plot = p2, width = 11, height = 7)

```

```{r}

#Calculating Correlation between Age and Educational Rank
cor(df$age, df$educ_rank_generation, use='complete.obs')
```

### Additive Models (With and Without Covariates)


```{r}

m1b <- glm(voted_c ~ educ_c +  educ_rank_generation  + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "White"))

m2b <- glm(voted_c ~  educ_c + educ_rank_generation + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Black"))

m3b <- glm(voted_c ~  educ_c + educ_rank_generation  + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Latino"))

m4b <- glm(voted_c ~ educ_c + educ_rank_generation + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Asian"))

#With controls
m1a <- glm(voted_c ~  educ_c + educ_rank_generation  + age + I(age^2) + sex + married + as.factor(year),
                           family = "binomial", 
                           data = subset(df, race_5 == "White"))

m2a <- glm(voted_c ~ educ_c + educ_rank_generation + age + I(age^2) + sex + married + as.factor(year),
                           family = "binomial", 
                           data = subset(df, race_5 == "Black"))

m3a <- glm(voted_c ~ educ_c + educ_rank_generation  + age + I(age^2) + sex + married + as.factor(year),
                             family = "binomial", 
                             data = subset(df, race_5 == "Latino"))

m4a <- glm(voted_c ~ educ_c + educ_rank_generation  + age + I(age^2) + sex + married + as.factor(year),
                             family = "binomial", 
                             data = subset(df, race_5 == "Asian"))

robust_se <- lapply(list(m1b, m2b, m3b, m4b, m1a, m2a, m3a, m4a), function(model) {
  round(sqrt(diag(vcovHC(model, type = "HC0"))), 3)
})
se = robust_se

stargazer(m1b, m2b, m3b, m4b, m1a, m2a, m3a, m4a,
          type = "html", 
          title = "Table 1: Effects of Educational Rank and Attaintment on Voter Turnout by Race/Ethnicity",
          covariate.labels = c("Education Attainment", "Education Rank", 
                               "Age", "Age Squared", "Female", "Married"),
          dep.var.labels = c("Voted (1 = Yes, 0 = No)"),
          omit = "year",
          add.lines = list(c("Year FE", rep("Yes", 8))),
          omit.stat = c("ll", "bic"),
          column.labels = c("White", "Black", "Latino", "Asian", "White", "Black", "Latino", "Asian"),
          coef = lapply(list(m1b, m2b, m3b, m4b, m1a, m2a, m3a, m4a), function(m) round(coef(m), 3)),
          se = robust_se,
          out = "models_new_add.htm")

        
```


#INTERACTIVE MODELS 

```{r}


m1d <- glm(voted_c ~ educ_c*educ_rank_generation  + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "White"))

m2d <- glm(voted_c ~  educ_c*educ_rank_generation + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Black"))

m3d <- glm(voted_c ~  educ_c*educ_rank_generation  + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Latino"))

m4d <- glm(voted_c ~ educ_c*educ_rank_generation + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Asian"))

# Model for White with interaction term
m1 <- glm(voted_c ~ educ_c*educ_rank_generation+ + age + I(age^2)  + sex +  married + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "White"))

# Model for Black race with interaction term
m2 <- glm(voted_c ~  educ_c*educ_rank_generation  + age + I(age^2) + sex + married + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Black"))

# Model for Latino race with interaction term
m3 <- glm(voted_c ~  educ_c*educ_rank_generation  + age + I(age^2) + sex + married + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Latino"))

#Model for Asian with interaction term
m4 <- glm(voted_c ~ educ_c*educ_rank_generation + age + I(age^2) + sex + married + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Asian"))

robust_se_interaction <- lapply(list(m1d, m2d, m3d, m4d, m1, m2, m3, m4), function(model) {
  round(sqrt(diag(vcovHC(model, type = "HC0"))), 3)
})

stargazer(m1d, m2d, m3d, m4d, m1, m2, m3, m4,
          type = "html", 
          title = "Table 2: Effects of Relative Education on Voter Turnout by Race/Ethnicity",
          dep.var.labels = c("Voted (1 = Yes, 0 = No)"),
          covariate.labels = c("Educ Attainment", "Educ Rank", 
                               "Age", "Age Squared", "Female", "Married", "Attainment X Rank", "Constant"),
          omit = "year",
          column.labels = c("White", "Black", "Latino", "Asian", "White", "Black", "Latino", "Asian"),
          add.lines = list(c("Year FE", rep("Yes", 8))),
          omit.stat = c("ll", "bic"),
          se = robust_se_interaction,
          out = "models_race_relative.htm")


```


#Marginal Effects Plots

```{r}

p1 <- interplot(m1,
                 var1 ='educ_c',
                 var2='educ_rank_generation',
                 hist = TRUE)  +
  ggtitle("White")+
  ylab("ME of Years of Education") +
  xlab("Educational Rank") 
p1

p2 <- interplot(m2,
                 var1 ='educ_c',
                 var2='educ_rank_generation',
                 hist = TRUE)  +
  ggtitle("Black")+
  ylab("ME of Years of Education") +
  xlab("Educational Rank") 
p2


p3 <- interplot(m3,
                 var1 ='educ_c',
                 var2='educ_rank_generation',
                 hist = TRUE)  +
  ggtitle("Latino")+
  ylab("ME of Years of Education") +
  xlab("Educational Rank") 

p3

p4 <- interplot(m4,
                 var1 ='educ_c',
                 var2='educ_rank_generation',
                 hist = TRUE)  +
  ggtitle("Asian")+
  ylab("ME of Years of Education") +
  xlab("Educational Rank") 
p4

library(patchwork)
final_layout <- (p1 + p2 + p3 + p4) 
  plot_layout(ncol = 2)
  
  final_layout
  
ggsave("interactions_all.png", plot = final_layout, width = 10, height = 6)


```

#Print Marginal Effects for Text Interpretation

```{r}

library(marginaleffects)
me1 <- predictions(m1, by=c("educ_c", "educ_rank_generation")) 
subset(me1, educ_rank_generation==0.4)

predictions(m1, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=12))
predictions(m1, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=18))
predictions(m2, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=12))
predictions(m2, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=18))
predictions(m3, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=12))
predictions(m3, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=18))
predictions(m4, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=12))
predictions(m4, newdata=datagrid(educ_rank_generation=c(0.4, 0.9), educ_c=18))

```



# Probability of Voting Figures


```{r}

predictions <- ggpredict(m1, terms = c("educ_c", "educ_rank_generation[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

# Plot the predictions with custom labels and legend title
pwhite <- plot(predictions) +
  labs(y = "Predicted Probability",
       x = "",
       color = "",
       title = "White") +
  scale_color_discrete(labels = c("Low educational rank", "High educational rank")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),  # Make y axis numbers larger and bold
    legend.position = "top"  # Move legend to the top of the plot
  ) +  
  scale_x_continuous(limits = c(10, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35, .85), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

# Display the plot
pwhite

```


```{r}
predictions <- ggpredict(m2, terms = c("educ_c", "educ_rank_generation[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

# Plot the predictions with custom labels and legend title
pblack <- plot(predictions) +
  labs(y = "",
       x = "",
       color = "",
       title = "Black") +
  scale_color_discrete(labels = c("Low educational rank", "High educational rank")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),
      legend.position = "none" ) + scale_x_continuous(limits = c(10, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35, .85), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

# Display the plot
pblack
```


```{r}

predictions <- ggpredict(m3, terms = c("educ_c", "educ_rank_generation[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

platino <- plot(predictions) +
  labs(y = "Predicted Probability",
       x = "Years of Education",
       color = "",
       title = "Latino") +
  scale_color_discrete(labels = c("Low educational rank", "High educational rank")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),
      legend.position = "none")  +  scale_x_continuous(limits = c(10, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35, .85), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

# Display the plot
platino

```


```{r}
# Asian

#Model for Asian with interaction term
m4 <- glm(voted_c ~ educ_c*educ_rank_generation + age + I(age^2) + sex + married + as.factor(year),
               family = "binomial", 
               data = subset(df, race_5 == "Asian"))

predictions <- ggpredict(m4, terms = c("educ_c", "educ_rank_generation[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

# Asian
pasian <- plot(predictions) +
  labs(y = "",
       x = "Years of Education",
       color = "",
       title = "Asian") +
  scale_color_discrete(labels = c("Low educational rank", "High educational rank")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),
      legend.position = "none")  +  scale_x_continuous(limits = c(10, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35, .85), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

# Display the plot
pasian
```

# Make figure arrangment
```{r}

final_layout <- (pwhite + pblack + platino + pasian) + 
  plot_layout(ncol =2 , nrow = 2)

final_layout

ggsave("pp_all.png", plot = final_layout, width = 12, height = 12)


```

# APPENDIX

# Section A: LR tests

```{r}
# Load necessary packages
library(sandwich)
library(knitr)
library(htmltools)

# Run likelihood ratio tests
lr_white <- anova(m1a, m1, test = "Chisq")
lr_black <- anova(m2a, m2, test = "Chisq")
lr_latino <- anova(m3a, m3, test = "Chisq")
lr_asian <- anova(m4a, m4, test = "Chisq")

# Create summary table
lr_table <- data.frame(
  Group = c("White", "Black", "Latino", "Asian"),
  Df = c(lr_white$Df[2], lr_black$Df[2], lr_latino$Df[2], lr_asian$Df[2]),
  Deviance = round(c(lr_white$Deviance[2], lr_black$Deviance[2], 
                     lr_latino$Deviance[2], lr_asian$Deviance[2]), 2),
  `p-value` = signif(c(lr_white$`Pr(>Chi)`[2], lr_black$`Pr(>Chi)`[2], 
                       lr_latino$`Pr(>Chi)`[2], lr_asian$`Pr(>Chi)`[2]), 3)
)

# Generate HTML table
html_output <- kable(lr_table, format = "html", table.attr = "style='width:60%;'", 
                     caption = "Likelihood Ratio Tests for Interaction Terms by Racial Group")

# Save as HTML file
save_html(HTML(html_output), file = "lr_test_results.html")

```



# Section B Rank (Within Racial Groups)

```{r}
m_white <- glm(voted_c ~ educ_rank_race* educ_c + age + I(age^2)  + sex + married + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "White"))

predictions <- ggpredict(m_white, terms = c("educ_c", "educ_rank_race[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

# Plot the predictions with custom labels and legend title
pwhite <- plot(predictions) +
  labs(y = "Predicted Probability",
       x = "",
       color = "",
       title = "White") +
  scale_color_discrete(labels = c("High education environment", "Low education environment")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),  # Make y axis numbers larger and bold
    legend.position = "top"  # Move legend to the top of the plot
  ) +  
  scale_x_continuous(limits = c(12, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35,1), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

# Display the plot
pwhite

m_black <- glm(voted_c ~ educ_rank_race* educ_c + age + I(age^2)  + sex + married + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "Black"))

predictions <- ggpredict(m_black, terms = c("educ_c", "educ_rank_race[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

# Plot the predictions with custom labels and legend title
pblack <- plot(predictions) +
  labs(y = "Predicted Probability",
       x = "",
       color = "",
       title = "Black") +
  scale_color_discrete(labels = c("High education environment", "Low education environment")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),
      legend.position = "none" ) + scale_x_continuous(limits = c(12, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35, .9), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

# Display the plot
pblack

m_latino <- glm(voted_c ~ educ_rank_race* educ_c + faminc_c + age + I(age^2)  + sex + married + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "Latino"))

predictions <- ggpredict(m_latino, terms = c("educ_c", "educ_rank_race[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

# Plot the predictions with custom labels and legend title
platino <- plot(predictions) +
  labs(y = "Predicted Probability",
       x = "Yers of Education",
       color = "",
       title = "Latino") +
  scale_color_discrete(labels = c("High education environment", "Low education environment")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),
      legend.position = "none")  +  scale_x_continuous(limits = c(12, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35, .9), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

# Display the plot
platino

m_asian <- glm(voted_c ~ educ_rank_race* educ_c  + age + I(age^2)  + sex + married + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "Asian"))


predictions <- ggpredict(m_asian, terms = c("educ_c", "educ_rank_race[.4,.9]"),
                         vcov_fun = "vcovHC", vcov_type = "HC0")

pasian <- plot(predictions) +
  labs(y = "Predicted Probability",
       x = "Years of Education",
       color = "",
       title = "Asian American") +
  scale_color_discrete(labels = c("High education environment", "Low education environment")) +
  theme_minimal() +
  theme(
    text = element_text(size = 16, face = "bold"),  # Adjust general font size and weight
    legend.title = element_text(size = 16, face = "bold"),  # Adjust legend title font size and weight
    plot.title = element_text(size = 18, face = "bold"),  # Make plot title larger and bold
    axis.text.x = element_text(size = 16, face = "bold"),  # Make x axis numbers larger and bold
    axis.text.y = element_text(size = 16, face = "bold"),
      legend.position = "top")  +  scale_x_continuous(limits = c(12, 18)) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0.35, .9), breaks = seq(0, 1, by = 0.1))  # Adjust y-axis to percentage scale

pasian
```






#Table B1: Effects of Relative Education within Racial Groups on Voter Turnout

```{r}
m_white <- glm(voted_c ~ educ_rank_race* educ_c  + age + I(age^2)  + sex + married + as.factor(year), family = "binomial", data = subset(df, race_5 == "White"))

m_black <- glm(voted_c ~ educ_rank_race* educ_c  + age + I(age^2)  + sex + married + as.factor(year), family = "binomial", data = subset(df, race_5 == "Black"))

m_latino <- glm(voted_c ~ educ_rank_race* educ_c  + age + I(age^2)  + sex + married + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "Latino"))

m_asian <- glm(voted_c ~ educ_rank_race* educ_c + age + I(age^2)  + sex + married + as.factor(year),
           family = "binomial", data = subset(df, race_5 == "Asian"))

stargazer(m_white, m_black, m_latino, m_asian)
stargazer(m_white, m_black, m_latino, m_asian,
          type = "html", 
          title = "Effects of Relative Education within Racial Groups on Voter Turnout",
          covariate.labels = c("Education Rank (Race)", "Education Attainment",
                               "Age", "Age Squared", "Female", "Married", 
                               "Education Rank (Race) x Education Attainment", "Constant"),
          omit = "year",  
            dep.var.labels = c("Voted (1 = Yes, 0 = No)"),
           column.labels = c("White", "Black", "Latino", "Asian"),
          add.lines = list(c("Year FE", "Yes", "Yes", "Yes", "Yes")),
          omit.stat = c("ll", "bic"),  out="models_racerank_relative.htm")
```

# Figure B1: Probability of Voting by Educational Rank (Race) and Attainment
```{r}
# Display the plot
final_layout_race <- (pwhite + pblack + platino + pasian) + 
  plot_layout(ncol = 2 , nrow = 2)
final_layout_race

ggsave("pp_all_rank_race.png", plot = final_layout_race, width = 14, height = 12)

```
















